#ifndef _MGGS_HERMITEINTERP_
#define _MGGS_HERMITEINTERP_

#include "InterpConditions.h"
#include "Polynomial.h"
#include <iostream>

using namespace std;

int frac[] = {1,2,6,24,120,720};

template <class T, int Num>
class HermiteInterp{
protected:
  InterpConditions<T,Num> PreConditions_;
  vector<vector<T> > tableOfDividedDiffs_;
  Polynomial<T> interPoly_;
public:
  HermiteInterp();
  HermiteInterp(const InterpConditions<T,Num>& IC);
  HermiteInterp(const string &_file);
  ~HermiteInterp();
  void setPreCondition(const InterpConditions<T,Num>& IC);
  void setPreCondition(const string &_file);
  const InterpConditions<T,Num>& getPreConditions() const;
  const vector<vector<T> >& getTableOfDividedDiffs() const;
  const Polynomial<T>& getInterPoly() const;
  //void showTableOfDividedDiffs() const;
  void Interpolation();
};


template <class T, int Num>
HermiteInterp<T,Num>::HermiteInterp(){};

template <class T, int Num>
HermiteInterp<T,Num>::HermiteInterp(const InterpConditions<T,Num>& IC){
  PreConditions_ = IC;
  this -> Interpolation();
}

template <class T, int Num>
HermiteInterp<T,Num>::HermiteInterp(const string &_file){
  this -> setPreCondition(_file);
  this -> Interpolation();
}

template <class T, int Num>
HermiteInterp<T,Num>::~HermiteInterp(){};


template <class T, int Num>
void HermiteInterp<T,Num>::setPreCondition(const InterpConditions<T,Num>& IC){
  PreConditions_ = IC;
}

template <class T, int Num>
void HermiteInterp<T,Num>::setPreCondition(const string &_file){
  PreConditions_.setTableByFile(_file);
}

template <class T, int Num>
const InterpConditions<T,Num>& HermiteInterp<T,Num>::getPreConditions() const{
 return PreConditions_;
}

template <class T, int Num>
const vector<vector<T> >& HermiteInterp<T,Num>::getTableOfDividedDiffs() const{
  return tableOfDividedDiffs_;
}

template <class T, int Num>
const Polynomial<T>& HermiteInterp<T,Num>::getInterPoly() const{
  return interPoly_;
}

// template <class T, int Num>
// void HermiteInterp<T,Num>::showTableOfDividedDiffs(const int N) const{
//   vector<vector<T> > adjustTable;
//   vector<vector<bool> > adjustBoolTable;
//   int M = tableOfDividedDiffs_.size();
//   vector<T> vec[M];
//   vector<bool> boolvec[M];
//   for (int i = 0 ; i < M ; i++){
//     vec[i].resize(M,0);
//     boolvec[i].resize(M,1);
//     if (i <= 1){
//          vec[i] = tableOfDividedDiffs_[i];
//          boolvec[i].resize(M,1);
//     }
//     else{
//          for (int j = 0 ; j < M ; j++){
//               if (j < (i - 1)){
//                    vec[i][j] = 0;
//                    boolvec[i][j] = 0;
//               }
//               else{
//                    vec[i][j] = tableOfDividedDiffs_[i][j - i + 1];
//                    boolvec[i][j] = 1;
//               }
//          }
//     }
//     adjustTable.push_back(vec[i]);
//     adjustBoolTable.push_back(boolvec[i]);
//   }
//   InterpConditions<T,5> tmpInterpCondition(adjustTable,adjustBoolTable);
//   cout << tmpInterpCondition;
// }

template <class T, int Num>
void HermiteInterp<T,Num>::Interpolation(){
  const InterpConditions<T,Num>& PreCondition = this->getPreConditions();
  const vector<vector<T> >& Table = PreCondition.getTable();
  const vector<vector<bool> >& BoolTable = PreCondition.getBoolTable();
  const int y = Table[0].size();
  const int x = Num + 1;
  vector<T> Tvec1;
  vector<T> Tvec2;
  for (int j = 0 ; j < y ; j++){
    for (int i = 1 ; i < x ; i++){
      if (BoolTable[i][j] == 1){
        Tvec1.push_back(Table[0][j]);
        Tvec2.push_back(Table[1][j]);  
      }
    }
  }
  if (!tableOfDividedDiffs_.empty())
    tableOfDividedDiffs_.clear();
  tableOfDividedDiffs_.push_back(Tvec1);
  tableOfDividedDiffs_.push_back(Tvec2);
  int M = Tvec1.size();
  vector<T> Tvec[M];
  for (int i = M - 1 ; i > 0 ; i--){
    Tvec[i].resize(i,0);
    int t = M - i;
    vector<T> Tvec2 = tableOfDividedDiffs_.back();
    for (int j = 0 ; j < i ; j++){
      T tmp = 0;
      if (Tvec1[j + t] != Tvec1[j])
        tmp = (Tvec2[j + 1] - Tvec2[j])*1.0/(Tvec1[j + t] - Tvec1[j]);
      else{
        auto it = find(Table[0].begin(),Table[0].end(),Tvec1[j]);
        int found = it - Table[0].begin();
        tmp = Table[t + 1][found]*1.0/frac[t - 1];
      }
      Tvec[i][j] = tmp;
    }
    tableOfDividedDiffs_.push_back(Tvec[i]);
  }
  vector<T> tmpvec;
  tmpvec.clear();
  tmpvec.push_back(1);
  Polynomial<T> tmpPolynomial(tmpvec);
  Polynomial<T> vecPoly[M];
  vector<T> vec4Poly[M];
  for (int i = 0 ; i < M ; i++){
    if (i > 0){
      vec4Poly[i].push_back(-Tvec1[i - 1]);
      vec4Poly[i].push_back(1);
      vecPoly[i].setCoe(vec4Poly[i]);
      tmpPolynomial = (tmpPolynomial * vecPoly[i]);
    }
    interPoly_ = interPoly_ + (tableOfDividedDiffs_[i+1][0] * tmpPolynomial);
  }
}


#else
//do nothing
#endif
